/*=====================================================================*
 *                   Copyright (C) 2011 Paul Mineiro                   *
 * All rights reserved.                                                *
 *                                                                     *
 * Redistribution and use in source and binary forms, with             *
 * or without modification, are permitted provided that the            *
 * following conditions are met:                                       *
 *                                                                     *
 *     * Redistributions of source code must retain the                *
 *     above copyright notice, this list of conditions and             *
 *     the following disclaimer.                                       *
 *                                                                     *
 *     * Redistributions in binary form must reproduce the             *
 *     above copyright notice, this list of conditions and             *
 *     the following disclaimer in the documentation and/or            *
 *     other materials provided with the distribution.                 *
 *                                                                     *
 *     * Neither the name of Paul Mineiro nor the names                *
 *     of other contributors may be used to endorse or promote         *
 *     products derived from this software without specific            *
 *     prior written permission.                                       *
 *                                                                     *
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND              *
 * CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES,         *
 * INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES               *
 * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE             *
 * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER               *
 * OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,                 *
 * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES            *
 * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE           *
 * GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR                *
 * BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF          *
 * LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT           *
 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY              *
 * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE             *
 * POSSIBILITY OF SUCH DAMAGE.                                         *
 *                                                                     *
 * Contact: Paul Mineiro <paul@mineiro.com>                            *
 *=====================================================================*/

#include <math.h>
#include <stdlib.h>
#include <string.h>

/* Ugh ... I should have used C++ */

#define logprating(_rating, _alpha, _logbeta, _n_labels)                \
  ({                                                                    \
    /*                                                                  
     * ok we want to compute the logarithm of
     * 
     * exp (exponents[rating])     exp (exponents[rating] - max)
     * ----------------------- = --------------------------------
     * sum_i exp (exponents[i])   sum_i exp (exponents[i] - max)
     *
     * the RHS is more numerically stable
     */                                                                 \
                                                                        \
    double _beta = exp (_logbeta);                                      \
    unsigned int _maxindex = 0;                                         \
                                                                        \
    for (unsigned int _logpri = 1; _logpri < _n_labels; ++_logpri)      \
      {                                                                 \
        if (-_alpha (_logpri) * _beta > -_alpha (_maxindex) * _beta)    \
          {                                                             \
            _maxindex = _logpri;                                        \
          }                                                             \
      }                                                                 \
                                                                        \
    double max = -_alpha (_maxindex) * _beta;                           \
    double denominator = 0;                                             \
                                                                        \
    for (unsigned int _logpri = 0; _logpri < _n_labels; ++_logpri)      \
      {                                                                 \
        denominator += exp (-_alpha (_logpri) * _beta - max);           \
      }                                                                 \
                                                                        \
    -_alpha (_rating) * _beta - max - log (denominator);                \
  })

#define dlogpratingdalpha(_dalpha, _rating, _alpha, _logbeta, _n_labels)      \
  ({                                                                          \
    double _beta = exp (_logbeta);                                            \
    unsigned int _maxindex = 0;                                               \
                                                                              \
    for (unsigned int _dlogprdai = 1; _dlogprdai < _n_labels; ++_dlogprdai)   \
      {                                                                       \
        if (-_alpha (_dlogprdai) * _beta > -_alpha (_maxindex) * _beta)       \
          {                                                                   \
            _maxindex = _dlogprdai;                                           \
          }                                                                   \
      }                                                                       \
                                                                              \
    double max = -_alpha (_maxindex) * _beta;                                 \
    double denominator = 0;                                                   \
    double ddenominator_dalpha[_n_labels];                                    \
                                                                              \
    memset (ddenominator_dalpha,                                              \
            0,                                                                \
            _n_labels * sizeof (ddenominator_dalpha[0]));                     \
                                                                              \
    for (unsigned int _dlogprdai = 0; _dlogprdai < _n_labels; ++_dlogprdai)   \
      {                                                                       \
        double d = -_alpha (_dlogprdai) * _beta - max;                        \
        double expd = exp (d);                                                \
        denominator += expd;                                                  \
        ddenominator_dalpha[_dlogprdai] -= expd * _beta;                      \
        ddenominator_dalpha[_maxindex] += expd * _beta;                       \
      }                                                                       \
                                                                              \
    for (unsigned int _dlogprdai = 0; _dlogprdai < _n_labels; ++_dlogprdai)   \
      {                                                                       \
        _dalpha[_dlogprdai] = -ddenominator_dalpha[_dlogprdai] / denominator; \
      }                                                                       \
                                                                              \
    _dalpha[_rating] -= _beta;                                                \
    _dalpha[_maxindex] += _beta;                                              \
  })

#define dlogpratingdlogbeta(_rating, _alpha, _logbeta, _n_labels)             \
  ({                                                                          \
    double _beta = exp (_logbeta);                                            \
    unsigned int _maxindex = 0;                                               \
                                                                              \
    for (unsigned int _dlogprdlbi = 1;                                        \
         _dlogprdlbi < _n_labels;                                             \
         ++_dlogprdlbi)                                                       \
      {                                                                       \
        if (-_alpha (_dlogprdlbi) * _beta > -_alpha (_maxindex) * _beta)      \
          {                                                                   \
            _maxindex = _dlogprdlbi;                                          \
          }                                                                   \
      }                                                                       \
                                                                              \
    double max = -_alpha (_maxindex) * _beta;                                 \
    double denominator = 0;                                                   \
    double ddenominator_dlogbeta = 0;                                         \
                                                                              \
    for (unsigned int _dlogprdlbi = 0;                                        \
         _dlogprdlbi < _n_labels;                                             \
         ++_dlogprdlbi)                                                       \
      {                                                                       \
        double d = -_alpha (_dlogprdlbi) * _beta - max;                       \
        double expd = exp (d);                                                \
        denominator += expd;                                                  \
        ddenominator_dlogbeta += expd * d;                                    \
      }                                                                       \
                                                                              \
    -_alpha (_rating) * _beta - max - ddenominator_dlogbeta / denominator;    \
  })
